home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / hlvector.lha / hl_vector / vali.cc < prev    next >
C/C++ Source or Header  |  1993-08-08  |  5KB  |  203 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. //
  3. //        Verification of the Aitken-Lagrange Interpolation
  4. //
  5.  
  6. #include "LinAlg.h"
  7. #include "math_num.h"
  8. #include <ostream.h>
  9.  
  10. /*
  11.  *------------------------------------------------------------------------
  12.  *        Set of the test functions to interpolate
  13.  */
  14.  
  15. static double (*F)(const double x);
  16.  
  17. static double f_exp(const double x)
  18. {
  19.   return exp(-x);
  20. }
  21.  
  22. static double f_sinexp(const double x)
  23. {
  24.   return sin(x) * exp(-x/10);
  25. }
  26.  
  27. /*
  28.  *------------------------------------------------------------------------
  29.  *    Set of tests for the interpolation over the uniform grid
  30.  */
  31.  
  32.             // Verify that interpolation exactly at a node
  33.             // really gives the function value at the node
  34. static void test_u_exact(const double s, const int ia, const int ib)
  35. {
  36.   cout << "\nCheck to see that interpolation at a node gives an exact"
  37.           " result\n";
  38.   cout << "\nUniform grid " << form("[%g:%g] with the step %g",s*ia,s*ib,s);
  39.  
  40.   Vector y(ia,ib);
  41.   Vector y_int(ia,ib);
  42.  
  43.   register int i;
  44.   for(i=y.q_lwb(); i<=y.q_upb(); i++)
  45.     y(i) = (*F)(s*i);
  46.  
  47.   for(i=y_int.q_lwb(); i<=y_int.q_upb(); i++)
  48.     y_int(i) = ali(s*i,s*ia,s,y);
  49.  
  50.   assert(y == y_int);
  51.  
  52.   cout << "\nDone\n";
  53. }
  54.  
  55.             // Check the precision for interpolation
  56.             // at arbitrary points
  57. static void test_u(const double s, const int ia, const int ib)
  58. {
  59.   cout << "\nCheck the precision of the interpolation at arbitrary points\n";
  60.   cout << "\nUniform grid " << form("[%g:%g] with the step %g",s*ia,s*ib,s);
  61.  
  62.   Vector y(ia,ib);
  63.   register int i;
  64.   for(i=y.q_lwb(); i<=y.q_upb(); i++)
  65.     y(i) = (*F)(s*i);
  66.  
  67.   cout << "\nPoint    Exact function value     Interpolated    Error\n";
  68.   register double q;
  69.   for(q=(ia-1)*s; q < (ib+2)*s; q += 1.1*s)
  70.   {
  71.     double yexact = (*F)(q);
  72.     double yint = ali(q,ia*s,s,y);
  73.     cout << form("%4.2g\t%12.6g\t\t%12.6g  %12.6g\n",q,yexact,yint,
  74.          abs(yexact-yint));
  75.   }
  76.  
  77.   cout << "\nDone\n";
  78. }
  79.  
  80. /*
  81.  *------------------------------------------------------------------------
  82.  *    Set of tests for the interpolation over the non-uniform grid
  83.  */
  84.  
  85.             // Verify that interpolation exactly at a node
  86.             // really gives the function value at the node
  87. static void test_n_exact(const Vector& x)
  88. {
  89.   cout << "\nCheck to see that interpolation at a node gives an exact"
  90.           " result\n";
  91.  
  92.   Vector y(x);
  93.   Vector y_int(x);
  94.  
  95.   register int i;
  96.   for(i=y.q_lwb(); i<=y.q_upb(); i++)
  97.     y(i) = (*F)(x(i));
  98.  
  99.   for(i=y_int.q_lwb(); i<=y_int.q_upb(); i++)
  100.     y_int(i) = ali(x(i),x,y);
  101.  
  102.   assert(y == y_int);
  103.  
  104.   cout << "\nDone\n";
  105. }
  106.  
  107.             // Check the precision for interpolation
  108.             // at arbitrary points
  109. static void test_n(const Vector& x,const double s, const int ia, const int ib)
  110. {
  111.   cout << "\nCheck the precision of the interpolation at arbitrary points\n";
  112.  
  113.   Vector y(x);
  114.   register int i;
  115.   for(i=y.q_lwb(); i<=y.q_upb(); i++)
  116.     y(i) = (*F)(x(i));
  117.  
  118.   cout << "\nPoint    Exact function value     Interpolated    Error\n";
  119.   register double q;
  120.   for(q=ia*s; q < ib*s; q += s)
  121.   {
  122.     double yexact = (*F)(q);
  123.     double yint = ali(q,x,y);
  124.     cout << form("%4.2g\t%12.6g\t\t%12.6g  %12.6g\n",q,yexact,yint,
  125.          abs(yexact-yint));
  126.   }
  127.  
  128.   cout << "\nDone\n";
  129. }
  130.  
  131.             // Check the precision for interpolation
  132.             // on example in the book
  133. static void test_n_book(void)
  134. {
  135.   cout << "\nCheck the precision of the interpolation at arbitrary points\n";
  136.   cout << "\nExample from Fig. 4.11 of the book"
  137.           "\n\tNumerical Methods and Software,"
  138.       "\n\tby D.Kahaner, C.Moler, and S.Nash - Prentice Hall, 1989\n";
  139.  
  140.   REAL xy[2][11] =
  141.   { {0,  2,  3,  5,  6,  8,  9,    11, 12, 14, 15},    // abscissae
  142.     {10, 10, 10, 10, 10, 10, 10.5, 15, 50, 60, 85}    // ordinates
  143.   };
  144.  
  145.   Vector x(0,sizeof(xy[0])/sizeof(xy[0][0])-1);
  146.   Vector y(x);
  147.   register int i;
  148.  
  149.   for(i=0; i<sizeof(xy[0])/sizeof(xy[0][0]); i++)
  150.     x(i) = xy[0][i], y(i) = xy[1][i];
  151.  
  152.   cout << "\n\t\tInterpolation nots";
  153.   cout << "\nx ";
  154.   for(i=x.q_lwb(); i<=x.q_upb(); i++)
  155.     cout << form("%6.2g ",x(i));
  156.   cout << "\ny ";
  157.   for(i=y.q_lwb(); i<=y.q_upb(); i++)
  158.     cout << form("%6.2g ",y(i));
  159.  
  160.   cout << "\n\nPoint    Interpolated value\n";
  161.   register double q;
  162.   for(q=0; q <= 16; q += 16/50.)
  163.     cout << form("%7.4g   %12.6g\n",q,ali(q,x,y));
  164.  
  165.   cout << "\nDone\n";
  166. }
  167.  
  168. /*
  169.  *------------------------------------------------------------------------
  170.  *                Root module
  171.  */
  172.  
  173. main()
  174. {
  175.   cout << "\n\n\t\tVerification of the Aitken-Lagrange interpolation\n";
  176.  
  177.   cout << "\n------------- Interpolation over the table with uniform grid";
  178.  
  179.   cout << "\nFunction to interpolate: exp(-x)\n";
  180.   F = f_exp;
  181.  
  182.   const double s = 0.1;
  183.   test_u_exact(s,1,10);
  184.   test_u(s,1,10);
  185.  
  186.   printf("\n------------- Interpolation over the table with non-uniform grid");
  187.  
  188.   printf("\nFunction to interpolate: sin(x) * exp(-x/10)\n");
  189.   F = f_sinexp;
  190.  
  191.   register int i;
  192.   Vector nodes(2,51);
  193.   for(i=nodes.q_lwb(); i<=11; i++)
  194.     nodes(i) = (i-1.)/10;
  195.   for(i=12; i<=nodes.q_upb(); i++)
  196.     nodes(i) = (i-6.)/5;
  197.   cout << "\nGrids 0.1 .. 1 have the mesh 0.1, and 1..9 have the mesh 0.2\n";
  198.  
  199.   test_n_exact(nodes);
  200.   test_n(nodes,0.27,0,15);
  201.   test_n_book();
  202. }
  203.